################################################################################
# Scandal Reaction
################################################################################


# Loading packages -------------------------------------------------------------
if(require(tidyverse) == F) {install.packages("tidyverse"); require(tidyverse)}
if(require(ggplot2) == F) {install.packages("ggplot2"); require(ggplot2)}
if(require(scales) == F) {install.packages("scales"); require(scales)}
if(require(here) == F) {install.packages("here"); require(here)}
if(require(lubridate) == F) {install.packages("lubridate"); require(lubridate)}
if(require(modelsummary) == F) {install.packages("modelsummary"); require(modelsummary)}
if(require(gt) == F) {install.packages("gt"); require(gt)}
if(require(did) == F) {install.packages("did"); require(did)}
if(require(ggeffects) == F) {install.packages("ggeffects"); require(ggeffects)}
if(require(scales) == F) {install.packages("scales"); require(scales)}
if(require(patchwork) == F) {install.packages("patchwork"); require(patchwork)}


### Set output directory -------------------------------------------------------
outfolder <- 'G:/Meu Drive/Papers/scandal_reaction/output'

################################################################################
##### Data
### loading data ---------------------------------------------------------------
load(file = file.path(outfolder, "scandal.reaction.testdata4.RData"))


################################################################################
### Subsample ------------------------------------------------------------------

# Function that selects municipalities where the party is electorally effective
# in the legislative branch across all periods before the first exogenous intervention
restriction <- function(party.label="PT"){
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(nep.is>0 & ANO_ELEICAO=='2004') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist()
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(nep.is>0 & ANO_ELEICAO=='2008') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist() %>% 
    intersect(restricao)
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(nep.is>0 & ANO_ELEICAO=='2012') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist() %>% 
    intersect(restricao)
  return(restricao)
}

# Function that selects municipalities where the party secures at least 0.1 of the
# councilor vote across all periods before the first exogenous intervention
restriction.seats <- function(level=.1, party.label="PT"){
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(vote.counc>level & ANO_ELEICAO=='2004') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist()
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(vote.counc>level & ANO_ELEICAO=='2008') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist() %>% 
    intersect(restricao)
  restricao <- mun.sample %>% filter(SIGLA_PARTIDO==party.label) %>% 
    filter(vote.counc>level & ANO_ELEICAO=='2012') %>% 
    select(COD_MUN_IBGE) %>% distinct() %>% unlist() %>% 
    intersect(restricao)
  return(restricao)
}


# Identify municipalities in which each party remained an effective party
# across all periods prior to the first exogenous treatment intervention:
mun.pt <- restriction(party.label ="PT")
mun.mdb <- restriction(party.label ="MDB")
mun.psdb <- restriction(party.label ="PSDB")
mun.prog <- restriction(party.label ="POGRESSISTAS")
mun.uni <- restriction(party.label ="UNIAO")
mun.ptb <- restriction(party.label ="PTB")
mun.pdt <- restriction(party.label ="PDT")
mun.pl <- restriction(party.label ="PL")
mun.psdb <- restriction(party.label ="PSDB")

# Select the sample of municipalities with effective parties:
# We must select municipalities where PT is effective, because PT (Workers' Party) 
# is the treatment group:
mun.sample.all <- mun.sample %>%
  filter(COD_MUN_IBGE %in% mun.pt) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.mdb) & SIGLA_PARTIDO=="MDB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.psdb) & SIGLA_PARTIDO=="PSDB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.prog) & SIGLA_PARTIDO=="POGRESSISTAS")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.uni) & SIGLA_PARTIDO=="UNIAO")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.ptb) & SIGLA_PARTIDO=="PTB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.pdt) & SIGLA_PARTIDO=="PDT")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.pl) & SIGLA_PARTIDO=="PL")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.psdb) & SIGLA_PARTIDO=="PSDB"))


# Identify municipalities in which each party held at least 0.1 of the seats
# in the legislature across all periods prior to the first exogenous treatment intervention
nivel=.1
mun.pt <- restriction.seats(level=nivel, party.label ="PT")
mun.mdb <- restriction.seats(level=nivel, party.label ="MDB")
mun.psdb <- restriction.seats(level=nivel, party.label ="PSDB")
mun.prog <- restriction.seats(level=nivel, party.label ="POGRESSISTAS")
mun.uni <- restriction.seats(level=nivel, party.label ="UNIAO")
mun.ptb <- restriction.seats(level=nivel, party.label ="PTB")
mun.pdt <- restriction.seats(level=nivel, party.label ="PDT")
mun.pl <- restriction.seats(level=nivel, party.label ="PL")
mun.psdb <- restriction.seats(level=nivel, party.label ="PSDB")

# Select the sample of municipalities with stable party presence:
# We must select municipalities where PT remains stable, because PT is the treatment group
mun.sample.seats <- mun.sample %>%
  filter(COD_MUN_IBGE %in% mun.pt) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.mdb) & SIGLA_PARTIDO=="MDB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.psdb) & SIGLA_PARTIDO=="PSDB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.prog) & SIGLA_PARTIDO=="POGRESSISTAS")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.uni) & SIGLA_PARTIDO=="UNIAO")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.ptb) & SIGLA_PARTIDO=="PTB")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.pdt) & SIGLA_PARTIDO=="PDT")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.pl) & SIGLA_PARTIDO=="PL")) %>% 
  filter(!(!(COD_MUN_IBGE %in% mun.psdb) & SIGLA_PARTIDO=="PSDB"))


################################################################################
##### Descriptive analysis

### Summary statistics for municipal data ---------------------------------------

# Select and prepare the appropriate data:
desc.sum <- mun.sample.all %>% 
  ungroup() %>% 
  # Create party dummies (one indicator per party):
  mutate(PT=ifelse(SIGLA_PARTIDO=='PT', 1, 0),
         MDB=ifelse(SIGLA_PARTIDO=='MDB', 1, 0),
         POGRESSISTAS=ifelse(SIGLA_PARTIDO=='POGRESSISTAS', 1, 0),
         PSB=ifelse(SIGLA_PARTIDO=='PSB', 1, 0),
         PSDB=ifelse(SIGLA_PARTIDO=='PSDB', 1, 0),
         PT=ifelse(SIGLA_PARTIDO=='PT', 1, 0),
         PTB=ifelse(SIGLA_PARTIDO=='PTB', 1, 0),
         PDT=ifelse(SIGLA_PARTIDO=='PDT', 1, 0),
         UNIAO=ifelse(SIGLA_PARTIDO=='UNIAO', 1, 0),
         PL=ifelse(SIGLA_PARTIDO=='PL', 1, 0)) %>% 
  mutate(post=ifelse(time>=2016, 1, 0)) %>%  # Post-treatment indicator: 1 if year ≥ 2016, else 0
  # Select treatment, time, outcomes, and covariates:
  select('post', "PT",
         "cand.mayor", "vote.mayor", 
         "cand.counc", "vote.counc", 
         # Covariates:
         "nep", "pib.def.log",
         "MDB", "POGRESSISTAS", "PSB", "PSDB", 
         "PTB", "PDT", "UNIAO", "PL")

# Human-readable variable labels:
new_names <- c("Post Treatment Period", "Work Party (Treated Group)", 
               "Mayoral Candidates", "Votes for Mayoral Candidates", 
               "Councilor Candidates", "Votes for Councilor Candidates", 
               "Number of Effective Paties", "Average Annual GDP (log)",
               "MDB", "POGRESSISTAS", "PSB", "PSDB", 
               "PTB", "PDT", "UNIAO", "PL")

# Build the final summary table:
desc.sum <- desc.sum %>%
  # Apply the new labels to the selected variables:
  rename_with(~set_names(new_names)) %>% 
  # Compute summary statistics with modelsummary::datasummary_skim:
  datasummary_skim(., fmt = 3, output = 'gt', histogram = F) %>% 
  # Add title and footnote:
  tab_header(title = md("Summary Statistics")) %>% 
  tab_footnote(footnote = "Notes: For each variable, the table displays the mean, 
  standard deviation, minimum, median, and maximum. 
               Each observation corresponds to a specific municipality.") %>% 
  # Improve table appearance:
  tab_options(heading.title.font.size=16) %>% 
  opt_vertical_padding(scale = 0) %>% 
  cols_hide(columns = c(2,3)) 

# Save to the output folder:
desc.sum %>% gtsave(filename = file.path(outfolder, "tab.summary.mun.tex"))

desc.sum

################################################################################
##### Tests

# Main results -----------------------------------------------------------------
# We restrict the sample to parties that were electorally effective and
# numerically relevant within the legislature across all pre-treatment periods.

# Effect on candidate nominations:
mod.cand.mayor <- att_gt(yname = "cand.mayor", gname = "treatment", 
                         idname = "id", tname = "time", 
                         data = mun.sample.all, base_period = "varying") 
mod.cand.counc <- att_gt(yname = "cand.counc", gname = "treatment", 
                         idname = "id", tname = "time",
                         data = mun.sample.all, base_period = "varying") 

# Effect on votes:
mod.vote.mayor <- att_gt(yname = "vote.mayor", gname = "treatment",
                         idname = "id", tname = "time", 
                         data = mun.sample.all, base_period = "varying") 
mod.vote.counc <- att_gt(yname = "vote.counc", gname = "treatment",
                         idname = "id", tname = "time", 
                         data = mun.sample.all, base_period = "varying") 

# Effect on votes including candidacies as endogenous covariates:
mod.vote.mayor.cand <- att_gt(yname = "vote.mayor", gname = "treatment",
                              idname = "id", tname = "time", xformla = ~cand.mayor,
                              data = mun.sample.all, base_period = "varying") 
mod.vote.counc.cand <- att_gt(yname = "vote.counc", gname = "treatment",
                              idname = "id", tname = "time", xformla = ~cand.counc,
                              data = mun.sample.all, base_period = "varying") 


# With time-variant covariates ------------------------------------------------------

# Select time-varying covariates:
covars <- c('nep', 'pib.def.log') %>% 
  paste(collapse = "+") %>% paste("~", .) %>% 
  as.formula()

# Effect on candidate nominations:
mod.cand.mayor.cov <- att_gt(yname = "cand.mayor", gname = "treatment", 
                             idname = "id", tname = "time", xformla = covars,
                             data = mun.sample.all, base_period = "varying") #xformla = covar.form, control_group=control
mod.cand.counc.cov <- att_gt(yname = "cand.counc", gname = "treatment", 
                             idname = "id", tname = "time", xformla = covars,
                             data = mun.sample.all, base_period = "varying") 

# Effect on votes:
mod.vote.mayor.cov <- att_gt(yname = "vote.mayor", gname = "treatment",
                             idname = "id", tname = "time", xformla = covars,
                             data = mun.sample.all, base_period = "varying") 
mod.vote.counc.cov <- att_gt(yname = "vote.counc", gname = "treatment",
                             idname = "id", tname = "time", xformla = covars,
                             data = mun.sample.all, base_period = "varying") 

# Effect on votes including candidacies as endogenous covariates:

# Add `cand.mayor` as an endogenous covariate:
covars <- c('nep', 'pib.def.log', 'cand.mayor') %>% 
  paste(collapse = "+") %>% paste("~", .) %>% as.formula()

mod.vote.mayor.cand.cov <- att_gt(yname = "vote.mayor", gname = "treatment",
                                  idname = "id", tname = "time", xformla = covars,
                                  data = mun.sample.all, base_period = "varying") 
mod.vote.counc.cand.cov <- att_gt(yname = "vote.counc", gname = "treatment",
                                  idname = "id", tname = "time", xformla = covars,
                                  data = mun.sample.all, base_period = "varying") 

# Restore the previous set of covariates—without the candidacy covariate:
covars <- c('nep', 'pib.def.log')


# Main results with seats' sample -----------------------------------------------
# We replicate the results using the sample that requires parties to have won
# at least 0.1 of legislative seats in all pre-treatment periods.

# Effect on candidate nominations:
mod.cand.mayor.seats <- att_gt(yname = "cand.mayor", gname = "treatment", 
                               idname = "id", tname = "time", 
                               data = mun.sample.seats, base_period = "varying") 
mod.cand.counc.seats <- att_gt(yname = "cand.counc", gname = "treatment", 
                               idname = "id", tname = "time",
                               data = mun.sample.seats, base_period = "varying") 

# Effect on votes:
mod.vote.mayor.seats <- att_gt(yname = "vote.mayor", gname = "treatment",
                               idname = "id", tname = "time", 
                               data = mun.sample.seats, base_period = "varying") 
mod.vote.counc.seats <- att_gt(yname = "vote.counc", gname = "treatment",
                               idname = "id", tname = "time", 
                               data = mun.sample.seats, base_period = "varying") 

# Effect on votes including candidacies as endogenous covariates:
mod.vote.mayor.cand.seats <- att_gt(yname = "vote.mayor", gname = "treatment",
                                    idname = "id", tname = "time", xformla = ~cand.mayor,
                                    data = mun.sample.seats, base_period = "varying") 
mod.vote.counc.cand.seats <- att_gt(yname = "vote.counc", gname = "treatment",
                                    idname = "id", tname = "time", xformla = ~cand.counc,
                                    data = mun.sample.seats, base_period = "varying") 

################################################################################
### Tables 
# Function table ---------------------------------------------------------------
# `get.did.tab` gets main results from a model list 
get.did.tab <- function(mod.list){
  # `mod.list` is a list with did models
  
  for(i in 1:length(mod.list)){
    # Getting agregated results:
    x <- mod.list[[i]] %>% aggte(type = 'simple')
    # option "simple" computes a weighted average of all group-time average treatment effects 
    # with weights proportional to group size
    
    out.tab <- data.frame(ATT=x$overall.att,
                          Std.Error=x$overall.se,
                          conf.low=x$overall.att-(1.96*x$overall.se),
                          conf.high = x$overall.att+(1.96*x$overall.se),
                          N=x$DIDparams$n,
                          conf.low.01=x$overall.att-(2.576*x$overall.se),
                          conf.high.01 = x$overall.att+(2.576*x$overall.se),
                          conf.low.001=x$overall.att-(3.291*x$overall.se),
                          conf.high.001 = x$overall.att+(3.291*x$overall.se)
    )
    # Note: Confidence intervals have gotten manually, as ATT - (1.96 * se), 
    # because values from `tidy` function do not correspond with confidence intervals
    # from `aggte` function (did package). To check and compare, run:
    # aggte(mod.gen.cand, type = 'simple')
    # tidy(mod.gen.cand)
    
    out.tab <- out.tab %>% 
      # Check statistical significance with at 95% confidence level:
      mutate(sig=ifelse(ATT>0 & conf.low>0 & conf.high>0, "*", ""),
             sig=ifelse(ATT<0 & conf.low<0 & conf.high<0, "*", sig),
             sig=ifelse(ATT>0 & conf.low.01>0 & conf.high.01>0, "**", sig),
             sig=ifelse(ATT<0 & conf.low.01<0 & conf.high.01<0, "**", sig),
             sig=ifelse(ATT>0 & conf.low.001>0 & conf.high.001>0, "***", sig),
             sig=ifelse(ATT<0 & conf.low.001<0 & conf.high.001<0, "***", sig))  %>% 
      # Rounding:
      mutate(ATT=round(ATT, 3), 
             ATT=paste0(ATT, sig),
             Std.Error=round(Std.Error, 3), 
             Std.Error=paste0("(", Std.Error, ")")) %>% 
      # Changing outcome variable lables:
      mutate(outvar=names(mod.list)[[i]]) %>% 
      select(outvar, ATT, Std.Error, N)
    
    if(i==1){did.out.tab <- out.tab}else{did.out.tab <- rbind(did.out.tab, out.tab)}
  }
  
  return(did.out.tab)
}

################################################################################
### Tables with time-variante covariates
### Mayoral Candidates -----------------------------------------------------------------
mod.list <- list(mod.cand.mayor.cov)
names(mod.list) <- c('mod.cand.mayor.cov')

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1)

# Replacing outcomes' lables:
names(tab) <- c("X", "Mayoral Candidates") #, "Councilor Candidates"

adlines <- data.frame(a=c("NEP", "GDP"),
                      b=c("Yes", "Yes"))
names(adlines) <- names(tab)
tab <- rbind(tab, adlines) %>% data.frame() 

# Getting Latex table:
tab.cand <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Impact of the Party-Linked Corruption Scandal on Candidate Nomination for Executive Offices') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on mayoral candidate nominations by parties, 
  using time-variant covariates and doubly robust estimation. 
  The coefficient is a weighted average of group-time treatment effects. 
  Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and 
  Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0) 

tab.cand %>% gtsave(filename = file.path(outfolder, "tab.cand.mayor.cov.tex"))

# Check the table:
tab.cand

### Councilos Candidates -----------------------------------------------------------------
mod.list <- list(mod.cand.counc.cov)
names(mod.list) <- c('mod.cand.counc.cov')

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1)

# Replacing outcomes' lables:
names(tab) <- c("X", "Councilor Candidates") #, "Councilor Candidates"

adlines <- data.frame(a=c("NEP", "GDP"),
                      b=c("Yes", "Yes"))
names(adlines) <- names(tab)
tab <- rbind(tab, adlines) %>% data.frame() 

# Getting Latex table:
tab.cand <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Impact of the Party-Linked Corruption Scandal on Candidate Nomination for Legislative Offices') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on councilor candidate nominations by parties, 
  using time-variant covariates and doubly robust estimation. 
  The coefficient is a weighted average of group-time treatment effects. 
  Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0) 

tab.cand %>% gtsave(filename = file.path(outfolder, "tab.cand.counc.cov.tex"))

# Check the table:
tab.cand

### Votes -------------------------------------------------------------
mod.list <- list(mod.vote.mayor.cov, mod.vote.mayor.cand.cov)
names(mod.list) <- c("mod.vote.mayor.cov", "mod.vote.mayor.cand.cov")

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1, V2)

# Replacing outcomes' lables:
names(tab) <- c("X", "Mayoral", " Mayoral")

adlines <- data.frame(a=c("Candidacy", "NEP", "GDP"),
                      b=c("No","Yes","Yes"),
                      c=c("Yes","Yes","Yes"))
names(adlines) <- names(tab)
tab <- rbind(tab, adlines) %>% data.frame() 

# Getting Latex table:
tab.votes <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Effect of Party-Linked Corruption Scandals on Party Vote Share') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on vote share of mayoral race by parties, 
  using time-variant covariates and doubly robust estimation. 
  The third column includes the proportion of candidacies as an endogenous covariate. 
  The coefficient is a weighted average of group-time treatment effects. 
  Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and 
  Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0)

tab.votes %>% gtsave(filename = file.path(outfolder, "tab.vote.cov.tex"))

# Check the table:
tab.votes


################################################################################
### Tables without time-variante covariates
### Mayoral Candidates -----------------------------------------------------------------
mod.list <- list(mod.cand.mayor)
names(mod.list) <- c('mod.cand.mayor')

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1)

# Replacing outcomes' lables:
names(tab) <- c("X", "Mayoral Candidates") #, "Councilor Candidates"

# Getting Latex table:
tab.cand <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Impact of the Party-Linked Corruption Scandal on Candidate Nomination for Executive Offices Without Time-Variante Covariates') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on mayoral candidate nominations by parties, 
  without time-variant covariates. The coefficient is a weighted average of 
  group-time treatment effects. Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and 
  Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0) 

tab.cand %>% gtsave(filename = file.path(outfolder, "tab.cand.mayor.tex"))

# Check the table:
tab.cand

### Councilos Candidates -----------------------------------------------------------------
mod.list <- list(mod.cand.counc)
names(mod.list) <- c('mod.cand.counc')

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1)

# Replacing outcomes' lables:
names(tab) <- c("X", "Councilor Candidates") #, "Councilor Candidates"

# Getting Latex table:
tab.cand <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Impact of the Party-Linked Corruption Scandal on Candidate Nomination for Legislative Offices Without Time-Variante Covariates') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on councilor candidate nominations by parties, 
  without time-variant covariates. The coefficient is a weighted average of 
  group-time treatment effects. Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0) 

tab.cand %>% gtsave(filename = file.path(outfolder, "tab.cand.counc.tex"))

# Check the table:
tab.cand

### Votes -------------------------------------------------------------
mod.list <- list(mod.vote.mayor, mod.vote.mayor.cand)
names(mod.list) <- c("mod.vote.mayor", "mod.vote.mayor.cand")

tab <- get.did.tab(mod.list)
tab <- tab %>% t() %>% as_tibble() 

# Selecting interest outcomes
tab <- tab[-1, ] %>% 
  mutate(V0=c('Corruption Scandal', "", "N")) %>% # Rename independent variables
  select(V0, V1, V2)

# Replacing outcomes' lables:
names(tab) <- c("X", "Mayoral", " Mayoral")

adlines <- data.frame(a=c("Candidacy"),
                      b=c("No"),
                      c=c("Yes"))
names(adlines) <- names(tab)
tab <- rbind(tab, adlines) %>% data.frame() 

# Getting Latex table:
tab.votes <- tab %>% 
  as.tbl() %>% 
  rename_with(~ str_replace_all(., "\\.|X", " ")) %>% 
  gt() %>%  
  tab_header(title = 'The Effect of Party-Linked Corruption Scandals on Party Vote Share Without Time-Variante Covariates') %>% 
  # Notes:
  tab_footnote(footnote = "Estimates show the average effect of 
  the party-linked corruption scandal on vote share of mayoral race by parties, 
  without time-variant covariates. The third column includes the proportion of 
  candidacies as an endogenous covariate. The coefficient is a weighted average of 
  group-time treatment effects. Standard errors are in parentheses. 
  These coefficients were estimated using the did R package by Callaway and 
  Sant’Anna (2021). 
  The number of observations appears at the table’s bottom.
  Significance levels are denoted as follows: * p < 0.05, ** p < 0.01, *** p < 0.001.") %>% 
  tab_options(heading.title.font.size = 16) %>% 
  opt_vertical_padding(scale = 0)

tab.votes %>% gtsave(filename = file.path(outfolder, "tab.vote.tex"))

# Check the table:
tab.votes



################################################################################
### Graphs 
# Graph function ---------------------------------------------------------------
get.did.graph <- function(model, ylabel, treat.title){
  
  model <- model %>% 
    aggte(., type='dynamic') %>% # Event-study aggregation: average ATT by event time (lengths of exposure relative to treatment)
    tidy() # Convert to a tidy tibble with estimate, conf.low, conf.high, event.time
  
  did.graph <- model %>% 
    ggplot(., aes(x = event.time, y = estimate)) +
    geom_point(aes(color = as.numeric(event.time)>=0), size = 2) +  # Points; color by pre (event.time<0) vs post (>=0)
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high, 
                      color = as.numeric(event.time)>=0), width = 0.2) + 
    scale_color_manual(values = c("TRUE" = "brown4", "FALSE" = "gray30"),
                       labels = c("Pre", "Post")) + # Manual legend for pre/post periods
    theme_minimal() +
    labs(x = "Years between elections", y = ylabel) + # Axis labels (y is passed in)
    ggtitle(treat.title) + # Panel title
    scale_x_continuous(breaks=seq(-12, 8, 4)) + # Event-time ticks (customizable)
    geom_hline(yintercept = 0, linetype = "dashed", color = "black") + # Zero line for visual reference
    geom_text(aes(label=paste0(round(model$estimate*100, 1), "%")), # Add % labels for point estimates
              position = position_dodge(1.4), size = 3.3, 
              vjust=-0.2, hjust=-0.05, colour = "black") +
    theme(axis.title = element_text(face = 'plain', colour="black", size=10),
          axis.ticks = element_blank(),
          axis.text.y = element_blank(),
          panel.grid.major.y = element_blank(),
          legend.position = 'bottom',
          legend.title = element_blank())
  
  if(length(model$estimate)==1){ # If only one time point, show a single legend label
    did.graph <- did.graph +
      scale_color_manual(values = c("TRUE" = "brown4", "FALSE" = "gray30"),
                         labels = c("Post"))
  }
  
  return(did.graph)
}

# Graphs -----------------------------------------------------------------------

g1.cand.mayor <- get.did.graph(mod.cand.mayor.cov, "Mayors", "")
g1.cand.counc <- get.did.graph(mod.cand.counc.cov, "Councilors", "")
g1.cand.mayor
ggsave(file.path(outfolder, 'g1.mayor.cand.png'), width = 7.5, height = 3.5)
g1.cand.counc
ggsave(file.path(outfolder, 'g1.counc.cand.png'), width = 7.5, height = 3.5)


g1.vote.mayor <- get.did.graph(mod.vote.mayor.cov, "Mayors", "Votes")
g1.vote.mayor.cand <- get.did.graph(mod.vote.mayor.cand.cov, "Mayors", "Votes with Candidacy Covariate")
(g1.vote.mayor + g1.vote.mayor.cand)
ggsave(file.path(outfolder, 'g1.votes.png'), width = 10, height = 7)


# Without covariates
g1.cand.mayor.ncov <- get.did.graph(mod.cand.mayor, "Mayors", "")
g1.cand.counc.ncov <- get.did.graph(mod.cand.counc, "Councilors", "")
g1.cand.mayor.ncov
ggsave(file.path(outfolder, 'g1.mayor.cand.ncov.png'), width = 7.5, height = 3.5)
g1.cand.counc.ncov
ggsave(file.path(outfolder, 'g1.counc.cand.ncov.png'), width = 7.5, height = 3.5)


g1.vote.mayor.ncov <- get.did.graph(mod.vote.mayor, "Mayors", "Votes")
g1.vote.mayor.cand.ncov <- get.did.graph(mod.vote.mayor.cand, "Mayors", "Votes without Candidacy Covariate")
(g1.vote.mayor.ncov + g1.vote.mayor.cand.ncov)
ggsave(file.path(outfolder, 'g1.votes.ncov.png'), width = 10, height = 7)


# Sample with parties with at least .10 of legislative seats
g1.cand.mayor.seats <- get.did.graph(mod.cand.mayor.seats, "Mayors", "")
g1.cand.counc.seats <- get.did.graph(mod.cand.counc.seats, "Councilors", "")
g1.cand.mayor.seats
ggsave(file.path(outfolder, 'g1.mayor.cand.seats.png'), width = 7.5, height = 3.5)
g1.cand.counc.seats
ggsave(file.path(outfolder, 'g1.counc.cand.seats.png'), width = 7.5, height = 3.5)


g1.vote.mayor.seats <- get.did.graph(mod.vote.mayor.seats, "Mayors", "Votes")
g1.vote.mayor.cand.seats <- get.did.graph(mod.vote.mayor.cand.seats, "Mayors", "Votes with Candidacy Covariate")
(g1.vote.mayor.seats + g1.vote.mayor.cand.seats)
ggsave(file.path(outfolder, 'g1.votes.seats.png'), width = 10, height = 7)


################################################################################
# Additional robustness checks
# One party at a time in control group -----------------------------------------

# Selecting parties from the control group:
control.list <- mun.sample.all$SIGLA_PARTIDO %>% unique() %>% 
  .[!grepl("^PT$", .)]

# Getting results for each party:
for(i in 1:length(control.list)){
  print(i)
  
  # Filter treatment and controle group:
    byparty.sample <- mun.sample.all %>%
      filter(SIGLA_PARTIDO==control.list[i] | SIGLA_PARTIDO=="PT")

    # Running models and getting:
    byparty.mayor <- att_gt(yname = "cand.mayor", gname = "treatment", 
                         idname = "id", tname = "time", #xformla = covars,
                         data = byparty.sample, base_period = "varying") %>% 
      # a weighted average of all group-time average treatment effects
      aggte(type = 'simple')

    if(i==1){byparty.mayor.g <- list()}
    byparty.mayor.g[[i]] <- att_gt(yname = "cand.mayor", gname = "treatment", 
                            idname = "id", tname = "time", #xformla = covars,
                            data = byparty.sample, base_period = "varying") %>% 
      get.did.graph(., "Votes", paste("PT X", control.list[i]))
    
    
    # Getting tables
    out.tab <- data.frame(ATT=byparty.mayor$overall.att,
                      Std.Error=byparty.mayor$overall.se,
                      conf.low=byparty.mayor$overall.att-(1.96*byparty.mayor$overall.se),
                      conf.high = byparty.mayor$overall.att+(1.96*byparty.mayor$overall.se),
                      N=byparty.mayor$DIDparams$n,
                      conf.low.01=byparty.mayor$overall.att-(2.576*byparty.mayor$overall.se),
                      conf.high.01 = byparty.mayor$overall.att+(2.576*byparty.mayor$overall.se),
                      conf.low.001=byparty.mayor$overall.att-(3.291*byparty.mayor$overall.se),
                      conf.high.001 = byparty.mayor$overall.att+(3.291*byparty.mayor$overall.se),
                      candidates="Mayors", contr.group=control.list[i])
    if(i==1){out.tab.mayor <- out.tab}else{out.tab.mayor <- rbind(out.tab.mayor, out.tab)}

}                      

# Binding tables:
byparty.tab <- rbind(out.tab.mayor) %>% 
  mutate(candidates="Mayoral Candidates") %>% 
  # For betther visualization:
  mutate(contr.group=gsub("PROGRESSISTAS", "PROGRES", contr.group),
         contr.group=paste("PT X", contr.group))

# Plot:
byparty.tab %>% 
  ggplot( ., aes(x = ATT, y = contr.group, xmin = conf.low, xmax = conf.high)) + 
  geom_errorbar(width = 0.002) +
  geom_point(size = 1.5) +
  facet_wrap( ~ candidates) +
  theme_minimal() +
  labs(x = 'ATT') +
  geom_text(aes(label=paste0(round(byparty.tab$ATT*100, 1), "%")), 
            position = position_dodge(1.4), 
            size = 3.7, vjust=-0.2, hjust=-0.05, colour = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") + 
  theme(axis.title.x = element_text(face = 'plain', colour="black", size=11),
        axis.title.y = element_blank(),
          axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = 'bottom',
        legend.title = element_blank())

ggsave(file.path(outfolder, 'by.party.cand.png'), width = 6, height = 4)

# Event Study Plots
(byparty.mayor.g[[1]] + byparty.mayor.g[[2]]) /
  (byparty.mayor.g[[3]] + byparty.mayor.g[[4]]) /
  (byparty.mayor.g[[5]] + byparty.mayor.g[[6]]) /
  (byparty.mayor.g[[7]] + byparty.mayor.g[[8]])

ggsave(file.path(outfolder, 'by.party.cand.trend.png'), width = 7, height = 10)


################################################################################
# By Number of Effective Parties -----------------------------------------------

# NEP average before the first treatment period
nep.av <- mun.sample.all %>%
  filter(ANO_ELEICAO=="2004"|ANO_ELEICAO=="2008"|ANO_ELEICAO=="2012") %>% 
  group_by(COD_MUN_IBGE) %>% summarise(nep.av=mean(nep)) %>% ungroup()

mun.sample.all <- mun.sample.all %>% 
  left_join(nep.av, by="COD_MUN_IBGE") %>% ungroup()

# Selecting parties from the control group:
nep.list <- mun.sample.all$nep.av %>% 
  quantile(na.rm = T, probs = seq(.05, 1, .05)) %>% 
  .[!grepl("^PT$", .)]

# Getting results for each party system:
for(i in 3:7){ 
  print(i)
  
  # Filter treatment and controle group:
  nep.sample <- mun.sample.all %>%
    filter(nep.av >= i & nep.av < (i+1))
  
  # Running models and getting:
  
  # For mayors
  covars <- c('nep', 'pib.def.log') %>% 
    paste(collapse = "+") %>% paste("~", .) %>% as.formula()
  byparty.mayor <- att_gt(yname = "cand.mayor", gname = "treatment", 
                          idname = "id", tname = "time",xformla = covars, 
                          data = nep.sample, base_period = "varying") %>% 
    # a weighted average of all group-time average treatment effects
    aggte(type = 'simple')
  
  if(i==3){byparty.mayor.g <- list()}
  byparty.mayor.g[[i]] <- att_gt(yname = "cand.mayor", gname = "treatment", 
         idname = "id", tname = "time",xformla = covars, 
         data = nep.sample, base_period = "varying") %>% 
    get.did.graph(., "Mayors", paste0(i,'≤NEP<',(i+1)))
  
  # Getting tables
  out.tab <- data.frame(ATT=byparty.mayor$overall.att,
                        Std.Error=byparty.mayor$overall.se,
                        conf.low=byparty.mayor$overall.att-(1.96*byparty.mayor$overall.se),
                        conf.high = byparty.mayor$overall.att+(1.96*byparty.mayor$overall.se),
                        N=byparty.mayor$DIDparams$n,
                        conf.low.01=byparty.mayor$overall.att-(2.576*byparty.mayor$overall.se),
                        conf.high.01 = byparty.mayor$overall.att+(2.576*byparty.mayor$overall.se),
                        conf.low.001=byparty.mayor$overall.att-(3.291*byparty.mayor$overall.se),
                        conf.high.001 = byparty.mayor$overall.att+(3.291*byparty.mayor$overall.se),
                        candidates="Mayors", nep=paste0(i,'≤NEP<',(i+1)))
  if(i==3){out.tab.mayor <- out.tab}else{out.tab.mayor <- rbind(out.tab.mayor, out.tab)}
  
}                      

# Binding tables:
out.tab.mayor <- out.tab.mayor %>% 
         mutate(nep=factor(nep, levels = rev(unique(out.tab.mayor$nep)))) 

# Plot
out.tab.mayor %>% 
  ggplot( ., aes(x = ATT, y = nep, xmin = conf.low, xmax = conf.high)) + 
  geom_errorbar(width = 0.002) +
  geom_point(size = 1.5) +
  theme_minimal() +
  labs(x = 'ATT') +
  geom_text(aes(label=paste0(round(out.tab.mayor$ATT*100, 1), "%")), 
            position = position_dodge(1.4), 
            size = 3.7, vjust=-0.2, hjust=-0.05, colour = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") + 
  theme(axis.title.x = element_text(face = 'plain', colour="black", size=11),
        axis.title.y = element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = 'bottom',
        legend.title = element_blank())

ggsave(file.path(outfolder, 'nep.cand.png'), width = 5, height = 4)

# Event Study Plots:
(byparty.mayor.g[[3]] + byparty.mayor.g[[4]]) /
(byparty.mayor.g[[5]] + byparty.mayor.g[[6]]) /
  byparty.mayor.g[[7]]

ggsave(file.path(outfolder, 'nep.cand.par.png'), width = 7, height = 8)

# End --------------------------------------------------------------------------
